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Abstract 

We consider the effect of intermolecular interactions on the optimal size- 
distribution of N hard spheres that occupy a fixed total volume. When 
we minimize the free-energy of this system, within the Percus-Yevick ap- 
proximation, we find that no solution exists beyond a quite low threshold ( 
r] ss 0.260). Monte Carlo simulations reveal that beyond this density, the 
size-distribution becomes bi-modal. Such distributions cannot be reproduced 
within the Percus-Yevick approximation. We present a theoretical argument 
that supports the occurrence of a non-monotonic size-distribution and em- 
phasizing the importance of finite size effects. 

I. INTRODUCTION 

Synthetic colloids are never perfectly monodisperse. Often, this polydispersity is a 
drawback — for instance, polydispersity is a problem in the preparation of high-quality col- 
loidal crystals, that are needed in photonic bandgap materials. However, occasionally, poly- 
dispersity is desirable, because it allows us to achieve material properties that cannot be 
realized with monodisperse colloids. For instance, monodisperse colloidal systems can fill 
at most 74.05 % of space in the crystalline phase (regular close packing) and some 63% in 
the liquid/glassy state (random close packing). In contrast, colloids with a properly chosen 
particle-size distribution can be made essentially space filling, both in the crystalline solid 
(Appolonian packing) and in the liquid. In practice, perfect space filling structures are never 
achieved because this requires an infinite number of (predominantly small) particles per unit 
volume. Here, we consider a somewhat simpler problem, namely the filling of a given volume 
V by a fixed number of particles N, that occupy a prescribed total volume fraction rj. We 
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assume that the particles are free to exchange volume. As we have fixed both the number 
and the total volume of the particles, the average volume per particle is fixed — it defines 
the natural length-scale in the model. Clearly, the Helmholtz free energy of the system will 
depend on the non-fixed particle-size distribution. The distribution, however, is restricted 
by the two constraints of fixed number of particles and fixed total volume. We define the 
optimal size distribution to be the one that minimizes the Helmholtz free energy under 
both constraints. In the Monte Carlo simulations (performed in the isothermal-isobaric en- 
semble) that we report in this paper, we study the density dependence of the particle-size 
distribution. We compare the simulation results for the size distribution with an analytical 
estimate that is obtained by solving the Percus-Yevick (PY) equation for an A-component 
hard-sphere mixture [|l[|. Not surprisingly, the PY equation works very well at low densi- 
ties. However, the theory breaks down at a surprisingly low density (77 w0.26). Of course, 
the fact that an approximate theory fails at a given density, does not imply that there is 
anything special going on in the system at that density. Yet, our simulations indicate that 
there is — the size distribution that was initially uni-modal, becomes bi-modal. We present 
a theoretical argument supporting this scenario. 

The remainder of this paper is organized as follows: in section |H| we describe the constant- 
pressure Monte Carlo simulations. The Percus-Yevick expression for the free energy of a 
system of polydisperse hard spheres with variable size distribution is discussed in section 
HI and |TV]. Section |V| is devoted to the derivation of analytical results turning useful when 
interpreting the simulation data. The mechanism behind the transition from uni-modal to 
bi-modal size distribution is also discussed. 

II. MONTE CARLO SIMULATIONS 

Monte Carlo simulations were performed in the isothermal-isobaric (constant NPT) 
ensemble ||. This means that the number of particles N, the pressure P and the temperature 
T are fixed. We attempt three distinct types of trial moves. We change the positions of the 
particles and allow the volume of simulation box to fluctuate, in order to equilibrate with 
respect to the applied pressure. Since we do not expect any crystalline order at low pressures, 
a cubic box shape is maintained. The third type of move is the one related to sampling the 
polydispersity of the system. To this end, we select two particles at random, between which 
we exchange an amount of volume drawn uniformly from the interval [— AV max , AV max ] (Fig. 
[I]). The maximum volume change AV max was chosen such that the acceptance of a volume 
exchange move is between 35 and 50 %. The relative frequency of the three moves is given 
by iV : 1 : N/2. The initial configurations are made by N monodisperse spheres on a simple 
cubic lattice. 

Simulations were performed for system sizes of iV=512 or 1000, and several different 
reduced pressures P* = (k B T)~ 1 P(cr 3 ), where a is the diameter of particles, and (•) denotes 
an average over particle-size distribution. As (a 3 ) is fixed, we choose (a 3 ) 1 / 3 to define 
the unit of length that we will use in the remainder of this paper. We use ksT as our 
unit of energy. All other units that we need, follow from these definitions. The equation 
of state [P* as a function of /?*(= p(c 3 ))] and the particle size distribution function were 
determined in the simulations. The results for the equation of state are shown in Fig. 
0. The simulation data have been collected in Table |. At low pressures the particle-size 
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distribution function is a single-peaked function with its maximum at v = (Figs. |3], |] 
and |5[). At higher pressures (typically, P* > 2.0), the particle- size distribution develops 
a second peak. Actually, this second peak is quite small (i.e. only a small fraction of all 
particles becomes "large"). However, these particles contribute appreciably to the total 
volume fraction (Fig. |6[). Depending on the pressure this contribution can get as large as 
75%. 

The formation of big particles in these MC simulations is a rather slow "dynamical" 
process. In order to speed up calculations, we did additional simulations, in which we 
started with a bi-disperse distribution, with one or several big particles containing 99% of 
the total volume occupied by the spheres, surrounded by a sea of small particles containing 
the remaining volume. In the 512 particle system, only one or two big particles remain 
for the lower pressures (P* = 2.5 and P* = 3.0). For higher pressures the number of big 
particles can stabilize at higher values as well. For the 1000 particle system the maximum 
number of big particles observed at the lower pressures is three. In addition, the size of 
these big particles is not the same. It is not clear whether this suggests a further possible 
fractionation or that it is a consequence of the slow equilibration and that one or more of 
the big particles are still shrinking. 

Below, we discuss these simulation results in the context of the relevant theoretical pre- 
dictions, but first let us stop to make certain considerations on the ideal entropy associated 
with this system. 



III. IDEAL ENTROPY OF A POLYDISPERSE SYSTEM 

Strictly speaking, the ideal entropy of a polydisperse system is infinite []!]. In a multi- 
component system such an entropy is exactly given by 

- NkB^wMA^pWi), (1) 

i 

where Wi is the molar fraction of species i and Aj its thermal wavelength. The usual path 
towards the entropy of a polydisperse system J^j (or for that matter, towards the entropy of 
continuous signals in Information Theory is to classify the species into "boxes" according 
to a certain property which distinguishes them (diameter, volume, molecular weight, etc). 
If x denotes such a property species i will denote the box having particles with x between 
iAx and [i + 1) Ax, for a given Ax which defines the boxes. If W(x) denotes the probability 
density of a particle having the value x for that property, then Wi = W(xi) Ax where Xi is a 
typical value of the zth box. Then Eq. ([]]) adopts the form 

- Nk B W(xi) Ax \n(A 3 tP W(xi)/Ax) . (2) 

i 

The entropy of the polydisperse system is obtained from Eq. (0) taking the limit A — > 0, 
and we can clearly see that, besides obtaining the usual expression 

Spoiy = —NK B [ dx W(x) ln[A(x) 3 W(x)} , (3) 



3 



there is also a divergent — In Ax, which is simply taken as a "constant" ignorable term. 

But expression (||) is not well defined. Suppose we simply change coordinates to label 
the species from x toy (say, from the diameter to the volume). Then the probability density 



in the new variable will be W(y) = W(x) 
labeling the entropy becomes 



It is straightforward to show that in the new 

dy 



S poly = -NK B J dyW(y) ln[A(x(y)) 3 W(y)] + Nk B J dyW(y) In 



dx 



(4) 



which is different of what we would have obtained had we began with the labeling y. 

This is a well known fact in Information Theory ||. In the study of fluid equilibria of 
quenched polydisperse systems this fact turns out to be irrelevant because the new term 
simply adds the same constant to both sides of the equilibrium equations [|J. However when 
studying annealed polydispersity this results tells us that the labeling is crucial and has to be 
dictated by the physical process underlying the polydispersity. In our case the Monte Carlo 
movements described in Section |TJ are a large scale description of a hypothetical microscopic 
system of tiny particles of exactly the same size distributed among N aggregates of a variable 
number of particles. The constant volume constraint would correspond to the conservation 
of the number of tiny particles, and the natural labeling of the aggregates would be the 
number of tiny particles which form it. As this number is proportional to the volume of the 
aggregate, in the continuum description is the volume, instead of the diameter, what turns 
out to be the natural labeling. 

Notice that we could have described another model in which the tiny particles aggregated 
into spherical surfaces. In that case it would be the total surface what would be preserved 
and the natural labeling of aggregates would be their respective surfaces. As we will discuss 
in the conclusions, the physics of this model would be completely different. 



IV. PERCUS-YEVICK THEORY 

The Percus-Yevick equation for an n-component hard-sphere mixture can be solved an- 
alytically, to yield the following equation of state ]l] : 

6 1-6 (i-£ 3 ) 2 (i-6) 3 ' {) 

The j-th moment of the particle-size distribution 6 is defined as 

i N ' 

where pi = Ni/V, the index i is used to denote the different particle species, and t>j is the 
volume of the z-th species. 

Equation (f|) is also valid for a continuous size distribution, in which case the sum in Eq. 
(§) is replaced by an integral. The corresponding expression for the chemical potential of a 
species with radius R is [|J 
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„• = In [ P A>W(v]] - ln(l - (,) 4- 4- ^ 4- + ^P'B? (7) 

(1-60 (1-6) (i-6) 3 

where A is the de-Broglie thermal wavelength, \/h 2 j (2irmkT), and W(v) is the probability 
density to find a particle with a volume around v = (4tt/3)R 3 . The pressure P* is given by 
Eq. (§. 

In an (NPT) description, the Gibbs free energy of the system fulfilling the constraints, 
must be at a minimum. The conservation of the number of particles and of the solid volume 
fraction, imply that W(v) must be of the form: 

W^)=exp j^a^j, (8) 

where 
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6 3£ ; 



" 2 = - 12( T^ + 2(T^- < W > 

The coefficients cto and «3 are determined by the constraints that the number of particles and 
the solid volume fraction are fixed. Note that all £j {i = 1, 2, 3) are positive. Moreover, £ 3 is 
equal to the volume fraction 77, and is therefore necessarily less than one. Hence, ot\ and a 2 
are always negative. The last coefficient, a^, should be negative or zero, because otherwise 
the particle-size distribution cannot be normalized. Since a%, a.2 and 03 are always negative 
the Percus-Yevick equation predicts that W(v) is a monotonically decreasing function of v. 
This implies that the size-distribution given by Eq. @ can never be bi-modal. 

Note that these conclusions also hold for the more accurate equation of state of Mansoori 
et al. 0. This equation adds an extra term to the pressure given in Eq. (|) depending on £ 3 . 
Thus Eq. ([Fp is the same with the new expression for P — which turns out to be irrelevant 
because the R 3 term is controlled by the Lagrange multiplier associated with the constraint 
on the total solid particle volume. In other respect, the analysis for the Mansoori equation- 
of-state is identical to that for PY. 

In practice, we solve Eq. (Q) numerically. To this end, we represent W(v) as a histogram. 
Initially, the value of W(v) in the different bins is assigned an arbitrary non-negative value, 
compatible with the constraint that W(v) is normalized and that (v) is fixed. We fix the 
density at the desired value. We determine the optimal W(v) using the following scheme: 
we select a bin (say i) at random and change the value of W(vi) by a random amount AW, 
distributed uniformly in the interval [— AW max , AW max \. We first check if the new value 
W(vi) is non-negative. If it is, we satisfy the constraints by scaling the width of all bins 
and the height of the function by two appropriately chosen factors. We then compute all 
moments £j, the pressure and the free energy, and we check if the Helmholtz free energy 
is smaller than the previous one. If it is, we accept the new value for W(vi), otherwise 
we reject it. We repeat the procedure until the free energy no longer decreases. We have 
verified that W(v) is indeed of the form given by Eqs. @ through fllDp . Figure | shows a 
comparison of the PY estimate for W(v), determined in this way, with the results of the full 
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Monte Carlo simulations. We find that «3 is a monotonically increasing function of density. 
A comparison between simulation and PY theory for the equation of state is shown in Fig. 
0. Note that, in this figure, the PY solution terminates at pressure P c ~ 1.34 (the cross in 
Fig. 0). This is the point where becomes zero. Beyond this point we can no longer find 
a solution for W(v) that is of the form given by Eq. (|8]). In Appendix [A], we consider the 
breakdown of the PY theory in more detail and obtain the packing fraction beyond which the 
PY approximation breaks down: 7] c = 0.260, the corresponding pressure being P* = 1.343. 
This breakdown of the PY equation at a relatively low density is surprising, as the PY 
equation works well up to quite high densities for fixed particle size distributions That 
it breaks down regardless the accuracy of the equation of state can be inferred from the fact 
that Mansoori et a/.'s equation of state undergoes exactly the same breakdown, though for 
slightly different values of r\ and P. Besides, from the analysis that we have carried out, it 
can be seen that a similar breakdown will appear for any other theory yielding an equation 
of state depending only on i = 0, 1, 2, 3. 



V. ANALYTICAL RESULTS 

In this section, we derive a theoretical bound for the pressure of the polydisperse system, 
providing the equation of state at the packing fraction where the size distribution becomes bi- 
modal. To this end, we work in the (NVT) ensemble and take advantage of the extremality of 
the Helmholtz free energy under the constraints of constant N and rj: the "grand potential" 

K = ^{W} - C J W(v) dv-Cx J vW{v) dv, (11) 

where Cq and C\ are Lagrange multipliers, has to be minimum for the optimal size distri- 
bution. In the above relation, the free energy functional T can be cast into the usual ideal 
and excess contributions 

T{W} = Nk B T J dv W(v) [In (A 3 pW(v)) - l] + F excess {W}. (12) 

We attempt the following change in the system: the volume of a given particle Vq is 
increased by an amount 5vq, before a rescaling of all volumes by a factor A (e.g. v — > Xv) 
such that the overall volume change vanishes. This imposes: 

The effect of the expansion of particle Vq on the size-distribution can be written: 

SW (v) = 1 [6{v-vq- 8v ) -S(v-v )}, (14) 
where 5(...) denotes the Dirac distribution. The scaling procedure affects W according to 

5W(v) = \ w (j) ~ W(v) (15) 

= W)^ + °^- (16) 



6 



The corresponding variation of the ideal contribution to T reads 



dv „.., ; 8W[v) 



with the functional derivative 

8W(v) 

We then get the entropic term 

ftFid = A;bT J dv ln[A 3 iy] 
/W>o) 



A; B T 



1 W(v ) + (v) 



5W(v) 



Nk B T ln[A 3 W»]. 



6(v-Vq- 5v ) - 5 (v - v ) + 



5Vn 



Svp d[vW] 
(v) dv 



(17) 
(18) 

(19) 
(20) 



where W is the derivative of W. 

The variation of the excess free energy reduces to the reversible work needed to perform 
the transformation, and is derived in appendix III: 



5v 



pk B T / dvW(v) g 



a + a 



1 



<ro J 



n 



(21) 



where g(ao/2 + cr/2) denotes the radial distribution function evaluated at contact between 
species of diameters o"o (having volume vq) and a (having volume v). When o"o (a), we 
can replace the density at the surface of particle Vq by that at a planar wall, and Eq. (|21~D 
becomes 



SW„ 



P 



n 



Sv n . 



In this limit Vq ^> (v) 

5v 5v 

For the optimal size distribution, 5F vanishes so that 



k B T 



f W'(v ) 1 



+ p 



W(V ) 1] 



2pk B T 
1 — r] 



P 



0. 



(22) 



(23) 



(24) 



Assuming W(v) to be a normalizable distribution, W'(v) has to be negative for large argu- 
ments, which sets the upper bound: 



P < 2 



pk B T 
l — i] 



or 



P* < 



12i] 



7T(1 -?7) 



(25) 



for the rescaled pressure. For low packing fractions (rj < i] c ) where the PY solution is 
available, the above inequality is fulfilled (Fig. 0). At the threshold 1] = i] c where the second 
polydispersity peak appears, W changes sign which means 
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The above relation is remarkably well obeyed within the PY approximation (see the data of 
section |IV] or figure 0: the PY expression crosses the line given by Eq. (j25l) exactly at rj c ). For 



rj > rjc, the upper bound is violated by the simulation results reported in Table | and figure |2|. 
However, the data suggest non negligible finite-size effects: increasing N shifts the pressure 
closer to the theoretical bound. Besides, starting from bi-disperse initial conditions (cf the 
procedure described in section 0), supposed to be closer to the expected polydispersity, has 
the same effect. According to expressions ( |2"B| ) and (|23), the violation of Eq. Q2"5D means 
that 51Z = bT < for 5v > 0, so that the biggest particle tends to expand. Its growth 
is however necessarily limited by the length L of the simulation box. This is supported 
by the observation that, even for the largest system investigated (N = 10 3 particles), the 
size of the biggest particle obtained is determined by L (o"bi gg est > L/3, irrespective of 
the packing fraction, see for example figure |6|). This suggests that system sizes that would 
presumably allow the system to reach thermodynamic equilibrium (and fulfill inequality (|25|) ) 
are numerically out of reach. Consequently, the question of the extensivity/intensivity of 
the number of large particles cannot be addressed by simulations; a theoretical investigation 
seems to require the detailed knowledge of the interfacial free energy between "large" and 
"small" species. At this stage, we cannot tell whether a true phase transition is associated 
with the occurrence of the second peak in the particle-size distribution. The available data 
certainly do not rule out this possibility. 

Finally, the integration of Eq. ( f24[) yields the tail of the optimal size-distribution: 



kiW(v) oc -v/(v) for v > (v). (27) 

Equation (||) obeys this relation, which cannot be tested against simulation results because 
of the lack of statistics for very large particles (not more than 5 in a typical run). 



VI. CONCLUSION 

At this stage, we can only speculate what will happen at larger N and/or larger densities. 
Conceivably, once the volume-fraction of the large particles exceeds a certain threshold, 
proliferation of still larger particles can occur, and so on, until eventually an "Appolonian" 
packing of the liquid is achieved. The theoretical analysis of this scenario is non-trivial, 
as the small particles now induce attractive depletion forces between the large particles. 
Unfortunately, the systems that we can conveniently study by simulation are too small to 
allow us to investigate this regime. 

We stress that the specific model we have chosen to study is somewhat arbitrary. For 
instance, rather than fixing the number of particles, one might have chosen to fix the total 
surface area of the particles. The latter constraint would be logical if one aims to model 
the size distribution of droplets covered with a fixed amount of surfactant. In addition, we 
assume that the surface free energy of the spheres is negligible. Again, this constraint can 
be removed. We hope that the rather surprising results of the present study will stimulate 
research into these related models. 
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APPENDIX A: 

Let us consider the range of densities where a solution of the Percus-Yevick equation is 
possible. As stated in section [TV], it is essential that be non-positive. Hence, the pressure 
at which 0:3 = defines the end point of the theory. To locate this point, consider the form 
of the solution at the point where a 3 = 0. Then the distribution function reduces to 



W(v) = exp I a i Ri 



i=0 



where R hereafter denotes a reduced radius measured in units of (<r 3 ) 1//3 . Using the two 
constraints for normalization and for the average volume of the particles, we can express the 
coefficients cto and a 2 in terms of ol\. If we combine Eqs. @ and fllPf) to eliminate £3, we 
obtain 

a 2 = 2«i-i - ~a\ (Al) 
£2 2 

(R) _i 2 



where the ratio of the moments £1/^2 does not depend explicitly on the density p or on «o, 
but it only depends on a\ itself. In the second line we have used the definition 

poo 

(if) = / R i W{v)AitR 2 dR. (A2) 
Jo 

But we have another relation between the moments of the particle-size distribution: partial 
integration of J R n W(v)dv yields 

POO 

[R n+1 W(v)]™ = / (2a 2 R n+2 + ai R n+1 + (n + l)R n ) W(v)dR = (A3) 
Jo 

for n > 0. This leads to the identities 



(flS> = 5 <fl2> -^ w - (A4) 

{&) = Z?1(R) - J-(R°). (A5) 
2a 2 2a 2 
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This allows us to write (R) and {R 2 ) as a function of a±, a 2 , (R°) and (-R 3 ); i.e. 

-3ai(iZ°) + 4cl<^ 3 ) 

^ = ^8^ ' (A6) 

= 12 W -2a 1 a aWj 
af — 8«2 



But = 1 and (R 3 ) = 1/8; hence 

(R) -3al + §a| 
W _ 12 - \ ai a 2 ' 

Substituting this expression in Eq. (|A1|) , we can eliminate (R) / (R 2 ) to obtain 



(A8) 



1 , — 3al + / » x 
a 2 + = A9 

2 12 - ^aia 2 

For what follows, it is convenient to introduce two new variables / and c 

<*i = -fc, (A10) 

"2 = -f; 



then 



( Rt ) = JT^( Rl )f=^ ( A11 ) 



We can use this relation to express / as function of c. We use the fact that the average 
volume per particle is fixed, to rewrite 

(R 3 ) _ 1 _ 1 (R 3 ) f=1 (c) 



{BP) 8 P(R°) f=1 (c 
1 

= 1 



- 3 Y(c) } (A12) 



where the second line defines the function Y(c). Hence 

/ = 2F 1/3 (c) (A13) 

As Y(c) can be expressed explicitly in terms of error functions, we now know / as a function 
of c. Equation ( |A10| ) allows us to express both a\ and a 2 as explicit functions of c. We 
can then use Eq. (|A9|) to determine c numerically. We find that this equation has a unique 
solution. Once the value of c has been determined, we know ai,a 2 and a (from the 
normalization condition). Eq. @ finally yields the packing fraction beyond which the PY 
approximation breaks down: i] c = 0.260198. The corresponding pressure is P* = 1.343442. 
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APPENDIX B: 



We first note that the reversible work done by an operator rescaling both particle volumes 
(vi — > Xvi, Vz) and container volume (V — > V = XV), is 

SW X = -Pidcai SV = -pk B T SV, where 5V = (X - 1) V. (Bl) 

In the transformation, the total volume V p of the particles changes according to 

f =f ^ IV, = iiV. (B2) 

Keeping the particle volumes fixed and going back to the original container volume (V — > 
V I X) requires the reversible work 

5W 2 — —P (V - V) = P SV. (B3) 

It is then straightforward to obtain the work associated with a rescaling of particle volumes 
at constant accessible volume V: 

5W v ^ Xv = 8W1 + 5W 2 = 8V P , (B4) 

rj 

valid for all polydispersities W{v). 

In the remainder, the shall derive the work needed to grow a particle Vq by an amount 
Svo (Svq = 7rcrQ(5cro/2). We assume the normalization J Wdv = 1 to hold. Consider species 
having volumes between v and v + Sv (diameters between a and a + 5a). They exert a 
pressure pkBTW(v)dv g(ao/2 + a/2) on particle vo, involving the radial distribution function 
at contact between species a and a. For the above pair (cr , a), the excluded volume sphere 
has diameter cr + a, and sweeps a volume 

SV 8weep = vr (a + a) 2 ^.= (l + ^ 5v (B5) 

during the growth of particle vq. Summing over all species v , the work performed by the 
operator takes the form 

a + a\ f a 



SW gTOWthV0 = pk B T 5v J dv W(v) g ) ( 1 + — ) • ( B 6) 

For the size modification considered in section M, the global volume change of the particles 
vanishes, such that 5V P + 5vo = 0. Summing the contributions arising from Eqs. ( B4|) and 
(|B6|), the results of Eq. (|21|) is recovered: 



^ = pkB T J dv W(v) g (2!±£) (l + £) * - . (B7) 



Note that a similar argument can be invoked to compute the reversible work needed 
to rescale all particle diameters, at constant V . After integration of Eq. (|B6|) over all 
W(vq) dv , we obtain: 
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8W v ^ Xv = pk B T j dvdv'W(v)W(v')g(^±^j {a + a'fa (B8) 



Inserting this result into Eq. (fFJ|) provides the equation of state for a polydisperse fluid of 
hard spheres 



P 



1+rjjdvdv' W{v) W{v') g (j^^j ip + °'Y 7^3 ■ ( B 9) 



pk B T 

For a monodisperse fluid, W(v) = S(v — vq) and we recover the well known relation 

JL- = l + 4 V g{ao). (BIO) 
pk B T 
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TABLES 



'/ 


p* 


0.0052(3) 


0.01 


0.244(1) 


1.00 


0.314(3) 


2.00 


0.374(2) 


2.50 


0.425(4) 


3.00 


0.484(2) 


4.00 


0.541(3) 


5.00 



TABLE I. Equation of state of polydisperse hard spheres, obtained from MC simulations with 
1000 particles. The estimated error in the last digit of the packing fraction rj is indicated in 
parentheses. 
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FIGURE CAPTIONS 



1. Schematic drawing of the Monte Carlo trial used to sample the polydispersity. Of two 
randomly chosen particles the volume of one is increased, while the other is decreased 
in volume by the same amount. 

2. The equation of state of the polydisperse system. The solid line is the PY prediction, 
which can not be extended beyond the cross. The circles correspond to a system 
initially prepared monodisperse, while the squares are final values in which we initially 
started with one big particle. The solid squares and the open symbols are from a 512 
and 1000 particle system respectively. The dashed line is the upper bound ( |26D for the 
pressure. 

3. The Monte Carlo results for the distribution of particle volumes, W(v), as a function 
of v, for several reduced pressures for a 1000 particle system. For P* > 2.0 the 
distribution develops a second peak (with statistical noise) at much larger volumes 
(note the change in scale on the right-hand side of the figure). Although there are 
only one or several of these big particles, they can contribute over 30% of the total 
volume of all particles. The diameters can get larger than a third of the length of the 
simulation box. 

4. Comparison of the numerical results (dashed line) for the particle-size distribution 
W(v), with the corresponding prediction of the PY theory (solid line). These results 
were obtained at a relatively low pressure (P* = 0.5), which corresponds to a volume 
fraction rj = 0.151. Note that at this density, simulation and theory are in quite good 
agreement. 

5. Snapshot of a typical configuration at reduced pressure P* = 1.0 and volume fraction 
rj = 0.244. 

6. Snapshot of a typical configuration at reduced pressure P* = 3.0 and volume fraction 
7] = 0.425. In this case there are two big particles. 
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FIG. 1. Zhang, Journal of Chemical Physics 
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FIG. 2. Zhang, Journal of Chemical Physics 
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FIG. 3. Zhang, Journal of Chemical Physics 
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FIG. 4. Zhang, Journal of Chemical Physics 
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FIG. 5. Zhang, Journal of Chemical Phy 
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FIG. 6. Zhang, Journal of Chemical Physics 
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